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ABSTRACT 

We develop a new nonlinear mean field dynamo theory that couples field growth 
to the time evolution of the magnetic helicity and the turbulent electromotive 
force, £. We show that the difference between kinetic and current helicities 
emerges naturally as the growth driver when the time derivative of £ is 
coupled into the theory. The solutions predict significant field growth in a 
kinematic phase and a saturation rate/strength that is magnetic Reynolds 
number dependent /independent in agreement with numerical simulations. The 
amplitude of early time oscillations provides a diagnostic for the closure. 

PACS codes: 95.30.Qd; 98.38.Am; 52.55.Ip, 52.30. Cv; 98.35.Eg; 96.60.Hv 

Introduction- Mean field dynamo (MFD) theory has been a useful framework for modeling 
the in situ origin of large-scale magnetic field growth in planets, stars, and galaxies [1-4], 
and has also been invoked to explain the sustenance of fields in fusion devices [5, 6]. 
However, whether the backreaction of the magnetic field itself prematurely quenches the 
MFD has been debated [7-24]. Recent progress has emerged from incorporating magnetic 
helicity evolution into the theory. 

To make explicit the problem to be solved, we first average the magnetic induction 
equation to obtain the basic MFD equation [1, 3]: 

d t B = Vx£ + Vx(VxB)-AV 2 B, (1) 

where B is the mean (large-scale) magnetic field in Alfven speed units, A = is the 
magnetic diffusivity in terms of the resistivity rj, V is the mean velocity which we set = 0, 
and £ = (v x b) is the turbulent electromotive force, a correlation between fluctuating 
velocity v and magnetic field b in Alfven units. Textbook treatments [1, 3] invoke 
£ = aB — /5VxB, where a and (3 are pseudoscalar and scalar correlations of turbulent 
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quantities respectively. In the kinematic theory [1] a = — (r c /3)(v • Vxv), where r c is a 
correlation time, and ([!]) is solved with a and (3 as input parameters. 

But the kinematic theory is incomplete. In a study of helical MHD turbulence, Ref. 18 
derived approximate evolution equations for the spectra of kinetic energy, magnetic energy, 
kinetic helicity, and magnetic helicity (= (A • Vx A), where B = Vx A). These calculations 
suggested that a ~ (r c /3)(b -Vxb — v-Vxv), the residual helicity, where () indicates 
spatial or ensemble average. This form has been employed in attempts to understand 
nonlinear dynamo quenching by coupling magnetic helicity conservation into the dynamo 
through the (b • Vxb) term [4,19-23] Although these studies [e.g. 4, 19, 20], wrote down 
an equation for the time evolution of a, they derived quenching formulae for a only for the 
steady- state. Only after a coupled nonlinear system of time-dependent large and small scale 
magnetic helicity equations were solved [15], was it apparent that a dynamical quenching 
model based on residual helicity reveals both a kinematic growth phase and an asymptotic 
resistively limited phase as seen in numerical experiments [12]. The dynamical approach 
has also been applied to dynamos with shear [16]. 

But even in these dynamical approaches, the £ was assumed to be proportional to the 
residual helicity. Here we show that the required residual helicity emerges not from £, but 
from d t £, and that including the d t £ equation in addition to the MFD and total magnetic 
helicity equations is essential for a complete MFD theory. We first derive d t £ and then 
derive the triplet of equations to be solved for the simple shear free helical dynamo whose 
solutions can be compared with existing numerical simulations. We discuss these solutions, 
physical implications, and the relation to previous work. 

Deriving dt£- A two-scale nonlinear quenching approach invoking the residual helicity 
in £ [15], captures the nonlinear dynamo saturation seen in simulations [12], but the 
derivation of the appropriate £ has been elusive. To couple (b • Vxb) of £ to the magnetic 
helicity conservation equation, the full small-scale field b must enter this correlation not a 
low order approximation. The essence of the puzzle [17] is that 

£{t) = ( v (t) x hit)) = f\v(t) x d f h(t'))dt' = - f\h(t) x d t rv(t'))dt', (2) 

Jo Jo 

assuming that 6(0) = and that t » so that (v(0) x h(t)) = 0. The last two terms are 
both exact expressions for £, however neither leads naturally the residual helicity entering 
a: upon using the induction equation for b, the second term on the right leads to a term 
oc /(v(£)Vxv(£))d£'. Using the Navier-Stokes equation for v, the last term contributes a 
term oc / (b(£)Vxb(£))<i£ / . One emerges with a choice rather than the dif ference between 
the two helicities. So how does the residual helicity emerge? 
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Rather than impose the form of the S, we solve for it dynamically using 

d t S= (d t vxb) + (vxd t b). (3) 
To proceed, we need equations for dtb and dtv. Assuming V ■ v = we have 

d t b = B ■ Vv — v ■ VB + V x ( v x b) — V x ( v x b) + AV 2 b, (4) 

and 

d t v g = P qi (B -Vbi + b- VBi - vVuj + (v- Vu<) + b • V6 ?; - (b • V^)) + uV\ + f g , 

(5) 

where f is a divergence- free forcing function uncorrelated with b, v is the viscosity, and 
P qi = (8gi — V _2 VqVj) is the projection operator that arises after taking the divergence 
of the incompressible Navier-Stokes equation to eliminate the total fluctuating pressure 
(magnetic + thermal). Using Reynolds rules [25] to interchange brackets with time and 
spatial derivatives, the 5th term of (^) and the 4th and 6th terms in the parentheses of (|5|) 
do not contribute when put into the averages so we ignore them. 

The contribution to d t S from the 3rd term in (H) can be derived by direct use of @ 
in configuration space. We assume isotropy of the resulting velocity and magnetic field 
correlations for terms linear in B. We also retain the triple correlations. The contribution 
to dtS from the 2nd term in (|]) also contributes terms linear in B, and triple correlations. 
Here the terms linear in B are best derived in Fourier space. For this, we follow the 
technique in the appendix of Ref. 22, which invokes the Fourier transform of the terms 
linear in B contributing to (d t v x b), supplemented by a linear expansion of the projection 
operator in k\ « &2, where k\ is the characteristic wavenumber of the bracketed or mean 
quantities and k% is the characteristic wavenumber of the fluctuating quantities b and v. 
Collecting all surviving terms, we then have for (Q) 

d t £ = |((b • Vxb) - (v ■ Vx v))B - |(v 2 )VxB + z/(V 2 v x b) + A(v x V 2 b) + T v + T M , 

(6) 

where T M = (v x Vx(v x b)) and Tj = ((e jqn (P qi (b ■ Vh - v • Vvi)b n ) are the triple 
correlations. Note that the 3rd, 4th, 6th and 8th terms in @ come from the (v x d t b) term 
of (§) and the 2nd 5th and 7th terms come from the (<9tv x b) term of (|j). 

We are primarily interested in the component of S parallel to B. For this we have 

d t £\\ = {{d t v x b) + (v x d t b)) ■ B/|B| + (v x b) • <9 t (B/|B|). (7) 
Substituting (|[) into (0) gives 

d t £\\ = aB 2 /|B| - f3B ■ VxB/|B| - Q£\\ (8) 
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where a = (l/3)((b • Vxb) — (v • Vxv)), ft = (l/3)(v 2 ), and ( accounts for microphysical 
dissipation terms, the last term of (0), and any additional contribution arising from 
T M + T v ^ 0. Note that a and ft appear similar to the usual a and ft dynamo coefficients in 
£, but they are fundamentally different because they are coefficients in d t £ (and thus have 
different units) and do not involve r c . Note also that if isotropy of like correlations were 
strongly violated, a and ft would be anisotropic tensors in analogy to tensor generalizations 
of a and ft [26, 27]. We have not considered that here. 

Following [15] we define large and small-scale magnetic helicities as Hf* = (A-B) vo i and 
H^ 1 = (a • h) vo i, where () vo i indicates a global spatial average. Then (b • Vxb) TOi = k\H^ 
and (B • VxB}„ ol = k\YL\ l . We define the small-scale kinetic helicity = (v ■ Vxv). 
We assume V = 0, and a force-free large-scale field for which the dtH^ 1 equation becomes 

2 1/2 1/2 

degenerate [15] with that of dtB . Then |B| = ki \Hi |. We can thus rewrite (pi) as 

d t £ ]{ = fcPKl^^Ff - H?)/3 - kT{H^/{\H^)x - (9) 



Dynamo Equations- We couple (|9|) to the equations for small and large-scale magnetic 
helicity evolution for a dynamo in which the kinetic energy is externally forced and V = 0. 
We interpret B, A and £ as the k\ (0 < k± < k 2 ) component of B, A and (v x b) of a closed 
system to facilitate comparison with simulations of Ref. 12. The total magnetic helicity, 
H M = (A ■ B) voh then satisfies [1] d t {A ■ B) vol = -2(E • B) voh where E = -<9 t A - V0, and 
(f) is the scalar potential. The large and small scale integrated magnetic helicity equations 
are then [12, 15, 28] 



d t (A ■ B) vol = 2(£ ■ B) wi - 2A(B • VxB) ro() (10) 

and 

8 t (a ■ b) vol = -2(£- B) vo i - 2A(b • Vxb) vol . (11) 
When B is force- free, the two-scale approximation allows us to write ( JTD| ) and ([TTj) as 

d t Hf J = 2£ ll k{ /2 \Hf I \ 1 / 2 - 2\k{H™ (12) 

and 

d t H^ = -2£\\k\ /2 \H™\ l l 2 - 2\klR™. (13) 

We need to solve fll2]), (0), and (0) after converting them into dimensionless form. We 
define the dimensionless quantities h\ = H^ 1 {k 2 /v 2 ) and h 2 = H^ 1 (k 2 /v 2 ), R m = v 2 /Xk 2 , 
Pr M = v/X, r = tv 2 k 2 , Q = ~£\\/v%, x = ft hi, C = C/v 2 k 2 = (1 + Pr M )/Rm, and use 
Hrf = —k 2 v\. For (0), (|T^) and (g) respectively, this gives 

d T ht = ^Qhl^ih/h) 1 / 2 - 2h 1 (k 1 /k 2 ) 2 /R m , (14) 
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8 T h 2 = 2Qh] /2 (k l /k 2 ) 1 ' 2 - 2h 2 /R mj (15) 

and 

d T Q = - (h/k 2 ) 1/2 h\ /2 (l + h 2 )/3 + (h/k 2 f 2 h{ /2 X - (Q- (16) 

Solutions- Since H 2 < 0, Hf 1 > and H^ 1 < for a growing solution. The solutions 
°f (li)> (15), and ( f[q) for two different i? m are shown in Figs. 1 & 2 over different time 



ranges for both ( = 2/R M and ( = 1 we use a oc ft in Fig. 1, but the resulting solutions 
are only weakly sensitive to the form of ft as shown in Fig. 2. In Figs.l&2 we also plotted 
the empirical fit formula to numerical simulations [12] (equation (54) of Ref 12) using our 
dimensionless parameters, to demonstrate the good agreement. 

In Figs. 1&2 we also compare the triplet solution of h± with the doublet solution [15] 
that results from solving ([14]) and ([15]) but with imposing £ = «B — /5VxB such that 
Q = Q d = — {k\jk 2 ) x l 2 }vl 2 {\ + h 2 )r c /3 + (ki/ k 2 ) 3 ^ 2 h]^ 2 XTc, where the correlation time r c is 
a free parameter taken to be ~ 1. Note that the present triplet solution does not involve r c 
in the dynamo coefficients a and ft. But a remarkable result emerges: Fig. 1 shows that the 
triplet solution matches the doublet solution at early times for £ = 1, which corresponds to 
a damping time ~ r c . This arises from a closure in which the triple correlations T M and 
T v lead to a damping with time constant ~ r c , and the damping suppresses the oscillations. 
Fig. 2 also shows that the ( = 2/R m and £ = 1 cases are indistinguishable at late times. 
We can also compare the kinematic regimes of the triplet and doublet solutions. The rise to 
the first peak of hi in Fig. la is independent of R m and there the two R m triplet solutions 
overlap. This is the kinematic regime. The end of the doublet kinematic regime occurs at 
h\ 1, as seen in Fig. 1. Again, the doublet and triplet match when ( ~ 1. In sum: The 
doublet solution emerges as the limit of the triplet solution when the triple correlations act 
as a damping term. This closure can be tested with future simulations. 

The maximum kinematic growth rate for hi is a function of ( because it occurs 
where Q is a minimum. If we ignore resistive terms so that h\ = —h 2 , and assume that 
C << 1 and x = 1/3, then Eq. fllBf ) implies that the maximum growth rate occurs 
when h\ ~ 1 — k\jk 2 . From Fig la, the minimum of Q during the first oscillation is 
~ —1/3 (found to be independent of ki/k 2 ). Setting d T hi ~ nh\, the maximum kinematic 
growth rate from (0) is then n ~ (2/3)(A; 1 /fc 2 ) 1 / 2 (l - h/k^- 1 ' 2 ~ 0.33, for k 2 = hk ± . 
However, when ( = 1, the minimum of Q from (|16|) occurs where Q = Qd- In this case, 
n ~ {2/2>){ki/k 2 ){l — ki/k 2 ) ~ 0.11 for k 2 = 5k±. This demonstrates that the kinematic 
growth rate for ( << 1 is ~ 3 times that for £ = 1, and why the triplet kinematic growth 
rate with £ = 1 matches that of the doublet. Both results are seen in Fig 1. 

Inspection of ([14]), ([IB]), and (|16|) reveals why there are oscillations for a positive seed 
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h\ and ( « 1 (and independent of whether Sit = 0) = or Sit = 0) ^ 0). As long as 
— 1 < h 2 < 0, Q grows more negative and hi and h 2 grow with mutually opposite signs. As 
h 2 passes through —1 from above, d T Q changes sign immediately but hi continues to grow 
positive, albeit more slowly, until Q changes sign. Then, d T hi changes sign and hi decreases. 
But d T h 2 changes sign when d T hi does, so when h 2 eventually passes back through —1 from 
below, d T Q reverses sign again, and eventually Q becomes negative and hi again grows. 
Large R m terms only weakly damp the oscillations. This describes what happens in Fig la. 
If instead, ( ~ 1, once d T Q is depleted by the growth of h 2 , the ( term of (|T6|) takes over 
and Q decays without oscillating. Then hi grows without oscillations (Fig lb). 

Discussion- Textbook kinematic MFD theories solve only the MFD equation itself [1]. 
Recent nonlinear approaches incorporating magnetic helicity evolution dynamically [15, 16] 
solve a doublet: the MFD equation (or the dtH^ 1 equation) and the total magnetic helicity 
evolution equation (the d t H^ + dfH^ equation). The present paper solves a triplet: the 
MFD equation, the total magnetic helicity evolution equation, and the d%S equation. Only 
the present approach shows how the difference between kinetic and current helicities (the 
residual helicity) emerges as the MFD driver in a time dependent theory. The residual 
helicity in turn couples to the total magnetic helicity evolution equation. The physical 
interpretation of the solutions for a closed system is that as the large scale helical field 
grows from MFD action, the small scale magnetic helicity grows of the opposite sign. At 
early times, kinematic growth is unimpeded, and the large scale field energy grows to 
B > [ki/k 2 )v 2 . Eventually, the small-scale magnetic helicity backreacts on the kinetic 
helicity, suppressing the growth rate to an R M dependent value. Ultimately B ~ (k 2 /ki)v 2 
at saturation. This picture also arises in the imposed S doublet approach [15]. The 
approaches agree in the asymptotic R m dependent growth phase, matching simulations 
[12]. However, oscillations are possible at early times only in the triplet approach. The 
amplitude of such oscillations serves as a direct diagnostic for the MHD closure scheme, 
which can be tested with future numerical experiments. For ( = 1, the agreement between 
the two approaches becomes exact for all times. This corresponds to the triple correlations 
contributing a simple damping term with characteristic time scale ~ r c . 
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Figure 1: (a) Plot for k± — 1, k 2 — 5 with ( = 2/R m . The top and middle oscillating 
curves are for hi with R m = 200 and R m = 1000, and the bottom oscillating curve 
is Q for R m = 1000. The thin lines are the doublet solutions for hi from Ref. 15 
which used an imposed S. Here % oc (1 + h 2 ). (b) Same as (a) but with £ ~ 1. 




Figure 2: (a) Same as Fig. la but for broader time range, (b) Same as (a) but for 
C — 1. For this time range, the doublet and triplet solutions are indistinguishable. 
The dotted curves are fits to simulations [12]. The lines slightly below each of the 
thick lines are for x oc 1/(1 + kihi/k 2 ) demonstrating the weak dependence on (3. 



